%% omega1 grid (equal spaced)

    % omega1 grid 
    temp_grid = linspace(-3*param.omega1_sd, 3*param.omega1_sd, param.n1+1);
    temp_lgrid = temp_grid(1:end-1);
    temp_rgrid = temp_grid(2:end);
    omega1_grid = (temp_lgrid'+temp_rgrid')/2;
    
    % omega1 density 
    temp_lcdf = normcdf(temp_lgrid, param.omega1_mean, param.omega1_sd);
    temp_rcdf = normcdf(temp_rgrid, param.omega1_mean, param.omega1_sd);
    omega1_density = temp_rcdf' - temp_lcdf'; 
    omega1_density = omega1_density/sum(omega1_density);

    
%% omega0 grid (equal weight)

    omega0_grid = icdf('normal', [0.1;0.3;0.5;0.7;0.9], param.omega0_mean, param.omega0_sd);
    omega0_density = ones(param.n0,1)/param.n0;
    
    
%% omega0, omega1, and joint density

    const.omega1 = repelem(omega1_grid, param.n0);
    const.omega0 = repmat(omega0_grid, param.n1, 1)+const.omega1*param.omega_rho; %introduce correlation    
    const.density_omega = repmat(omega0_density, param.n1, 1).*repelem(omega1_density, param.n0);

    
    
clearvars temp_* omega1_grid omega1_density omega0_grid omega0_density
    
    